Working with SWOT Level 2 Water Mask Raster Image Data Product:

Working with SWOT Level 2 Water Mask Raster Image Data Product:

In AWS Cloud Version

Author: Nicholas Tarpinian, PO.DAAC

Summary & Learning Objectives

Notebook showcasing how to work with SWOT Raster Image data product version 2.0 on a local machine

  • Utilizing the earthaccess Python package. For more information visit: https://nsidc.github.io/earthaccess/
  • Option to query the new dataset based on users choice; choosing between two resolutions either by ‘100m’ or ‘250m’.
  • Visualizing the dataset by choosing one of its various variables.
  • Visualizing multiple raster images on a single map.
  • Stacking multiple raster images and creating a time dimenstion to analyze over time.

Requirements

1. Compute environment

This tutorial is written to run in the following environment: - AWS instance running in us-west-2: NASA Earthdata Cloud data in S3 can be directly accessed via and s3fs session; this access is limited to requests made within the US West (Oregon) (code: us-west-2) AWS region.

2. Earthdata Login

An Earthdata Login account is required to access data, as well as discover restricted data, from the NASA Earthdata system. Thus, to access NASA data, you need Earthdata Login. Please visit https://urs.earthdata.nasa.gov to register and manage your Earthdata Login account. This account is free to create and only takes a moment to set up.

Import libraries

%matplotlib inline

import os
import io
import s3fs
import xarray as xr
import numpy as np
from datetime import datetime
from os.path import isfile, basename, abspath
from pathlib import Path
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.animation import FuncAnimation, PillowWriter
from IPython.display import display, Image, HTML
import hvplot
import hvplot.xarray 
import holoviews as hv
from holoviews import opts
import earthaccess
from earthaccess import Auth, DataCollections, DataGranules, Store

hv.extension('bokeh')

Authentication with earthaccess

In this notebook, we will be calling the authentication in the below cell.

auth = earthaccess.login(strategy="interactive", persist=True)

Search using earthaccess for SWOT Raster products

Each dataset has it’s own unique collection concept ID. For the SWOT_L2_HR_Raster_2.0 dataset, we can find the collection ID here.

For this tutorial, we are looking at the Lake Mead Reservoir.

We used bbox finder to get the exact coordinates for our area of interest.

raster_results = earthaccess.search_data(
    short_name = 'SWOT_L2_HR_RASTER_2.0',
    bounding_box=(-115.112686,35.740939,-114.224167,36.937819),
    temporal =('2024-01-01 00:00:00', '2024-01-05 23:59:59'),
    granule_name = '*_100m_*',
    count =200
)
Granules found: 2
fs_s3 = earthaccess.get_s3fs_session(results=raster_results)

Visualizing the Dataset

Let’s now visualize an individual layer for a single file path using Xarray to read the NetCDF image.

raster_link = earthaccess.results.DataGranule.data_links(raster_results[1], access='direct')[0]
s3_file_obj5 = fs_s3.open(raster_link, mode='rb')
ds_raster = xr.open_dataset(s3_file_obj5, engine='h5netcdf')
ds_raster
<xarray.Dataset>
Dimensions:                  (x: 1505, y: 1506)
Coordinates:
  * x                        (x) float64 5.931e+05 5.932e+05 ... 7.435e+05
  * y                        (y) float64 4.026e+06 4.026e+06 ... 4.176e+06
Data variables: (12/39)
    crs                      object ...
    longitude                (y, x) float64 ...
    latitude                 (y, x) float64 ...
    wse                      (y, x) float32 ...
    wse_qual                 (y, x) float32 ...
    wse_qual_bitwise         (y, x) float64 ...
    ...                       ...
    load_tide_fes            (y, x) float32 ...
    load_tide_got            (y, x) float32 ...
    pole_tide                (y, x) float32 ...
    model_dry_tropo_cor      (y, x) float32 ...
    model_wet_tropo_cor      (y, x) float32 ...
    iono_cor_gim_ka          (y, x) float32 ...
Attributes: (12/49)
    Conventions:                   CF-1.7
    title:                         Level 2 KaRIn High Rate Raster Data Product
    source:                        Ka-band radar interferometer
    history:                       2024-01-05T10:46:43Z : Creation
    platform:                      SWOT
    references:                    V1.1.1
    ...                            ...
    x_min:                         593100.0
    x_max:                         743500.0
    y_min:                         4025600.0
    y_max:                         4176100.0
    institution:                   CNES
    product_version:               01
ds_raster.wse.hvplot.image(y='y', x='x').opts(width=1000, height=600, colorbar=True)

Visualizing Multiple Images

Let’s now visualize the multiple raster images that were searched and explore the data.

Utilizing xarray.open_mfdataset which supports the opening of multiple files.

ds = xr.open_mfdataset(earthaccess.open(raster_results), engine='h5netcdf',combine='nested', concat_dim='x')
ds
Opening 2 granules, approx size: 0.06 GB
using endpoint: https://archive.swot.podaac.earthdata.nasa.gov/s3credentials
<xarray.Dataset>
Dimensions:                  (y: 2760, x: 3009)
Coordinates:
  * y                        (y) float64 3.9e+06 3.9e+06 ... 4.176e+06 4.176e+06
  * x                        (x) float64 5.683e+05 5.684e+05 ... 7.435e+05
Data variables: (12/39)
    crs                      (x) object b'1' b'1' b'1' b'1' ... b'1' b'1' b'1'
    longitude                (y, x) float64 dask.array<chunksize=(502, 502), meta=np.ndarray>
    latitude                 (y, x) float64 dask.array<chunksize=(502, 502), meta=np.ndarray>
    wse                      (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    wse_qual                 (y, x) float32 dask.array<chunksize=(2760, 1504), meta=np.ndarray>
    wse_qual_bitwise         (y, x) float64 dask.array<chunksize=(752, 752), meta=np.ndarray>
    ...                       ...
    load_tide_fes            (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    load_tide_got            (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    pole_tide                (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    model_dry_tropo_cor      (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    model_wet_tropo_cor      (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    iono_cor_gim_ka          (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
Attributes: (12/49)
    Conventions:                   CF-1.7
    title:                         Level 2 KaRIn High Rate Raster Data Product
    source:                        Ka-band radar interferometer
    history:                       2024-01-05T11:01:19Z : Creation
    platform:                      SWOT
    references:                    V1.1.1
    ...                            ...
    x_min:                         568300.0
    x_max:                         718600.0
    y_min:                         3900200.0
    y_max:                         4050500.0
    institution:                   CNES
    product_version:               01
    • y
      PandasIndex
      PandasIndex(Float64Index([3900200.0, 3900300.0, 3900400.0, 3900500.0, 3900600.0, 3900700.0,
              3900800.0, 3900900.0, 3901000.0, 3901100.0,
              ...
              4175200.0, 4175300.0, 4175400.0, 4175500.0, 4175600.0, 4175700.0,
              4175800.0, 4175900.0, 4176000.0, 4176100.0],
             dtype='float64', name='y', length=2760))
    • x
      PandasIndex
      PandasIndex(Float64Index([568300.0, 568400.0, 568500.0, 568600.0, 568700.0, 568800.0,
              568900.0, 569000.0, 569100.0, 569200.0,
              ...
              742600.0, 742700.0, 742800.0, 742900.0, 743000.0, 743100.0,
              743200.0, 743300.0, 743400.0, 743500.0],
             dtype='float64', name='x', length=3009))
  • Conventions :
    CF-1.7
    title :
    Level 2 KaRIn High Rate Raster Data Product
    source :
    Ka-band radar interferometer
    history :
    2024-01-05T11:01:19Z : Creation
    platform :
    SWOT
    references :
    V1.1.1
    reference_document :
    JPL D-56416 - Revision B - August 17, 2023
    contact :
    alexander.t.corben[at]jpl.nasa.gov
    cycle_number :
    8
    pass_number :
    511
    scene_number :
    109
    tile_numbers :
    [217 218 219 217 218 219]
    tile_names :
    511_217L, 511_218L, 511_219L, 511_217R, 511_218R, 511_219R
    tile_polarizations :
    V, V, V, H, H, H
    coordinate_reference_system :
    Universal Transverse Mercator
    resolution :
    100.0
    short_name :
    L2_HR_Raster
    descriptor_string :
    100m_UTM11S_N_x_x_x
    crid :
    PIC0
    pge_name :
    PGE_L2_HR_RASTER
    pge_version :
    5.0.4
    time_granule_start :
    2024-01-01T12:43:28.040906Z
    time_granule_end :
    2024-01-01T12:43:47.497154Z
    time_coverage_start :
    2024-01-01T12:43:28.041873Z
    time_coverage_end :
    2024-01-01T12:43:46.958982Z
    geospatial_lon_min :
    -116.24753568370588
    geospatial_lon_max :
    -114.56382182347716
    geospatial_lat_min :
    35.225727504637376
    geospatial_lat_max :
    36.59518831179102
    left_first_longitude :
    -116.24753568370588
    left_first_latitude :
    35.46460890403489
    left_last_longitude :
    -115.95968454013911
    left_last_latitude :
    36.59518831179102
    right_first_longitude :
    -114.86994488403649
    right_first_latitude :
    35.225727504637376
    right_last_longitude :
    -114.56382182347716
    right_last_latitude :
    36.35099280753312
    xref_l2_hr_pixc_files :
    SWOT_L2_HR_PIXC_008_511_217L_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_218L_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_219L_20240101T124346_20240101T124357_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_217R_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_218R_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_219R_20240101T124346_20240101T124357_PIC0_01.nc
    xref_l2_hr_pixcvec_files :
    SWOT_L2_HR_PIXCVec_008_511_217L_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_218L_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_219L_20240101T124346_20240101T124357_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_217R_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_218R_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_219R_20240101T124346_20240101T124357_PIC0_01.nc
    xref_param_l2_hr_raster_file :
    SWOT_Param_L2_HR_Raster_20000101T000000_21000101T000000_20230817T100000_v302.rdf
    xref_reforbittrack_files :
    SWOT_RefOrbitTrackTileBoundary_Nom_20000101T000000_21000101T000000_20200617T193054_v101.txt, SWOT_RefOrbitTrack125mPass1_Nom_20000101T000000_21000101T000000_20200617T193054_v101.txt, SWOT_RefOrbitTrack125mPass2_Nom_20000101T000000_21000101T000000_20200617T193054_v101.txt
    utm_zone_num :
    11
    mgrs_latitude_band :
    S
    x_min :
    568300.0
    x_max :
    718600.0
    y_min :
    3900200.0
    y_max :
    4050500.0
    institution :
    CNES
    product_version :
    01
    • # Lisiting all available variables to select and visualize.
variables = list(ds.data_vars)

print("Variables within the NetCDF file:")
for variable in variables:
    print(variable)
    Variables within the NetCDF file:
crs
longitude
latitude
wse
wse_qual
wse_qual_bitwise
wse_uncert
water_area
water_area_qual
water_area_qual_bitwise
water_area_uncert
water_frac
water_frac_uncert
sig0
sig0_qual
sig0_qual_bitwise
sig0_uncert
inc
cross_track
illumination_time
illumination_time_tai
n_wse_pix
n_water_area_pix
n_sig0_pix
n_other_pix
dark_frac
ice_clim_flag
ice_dyn_flag
layover_impact
sig0_cor_atmos_model
height_cor_xover
geoid
solid_earth_tide
load_tide_fes
load_tide_got
pole_tide
model_dry_tropo_cor
model_wet_tropo_cor
iono_cor_gim_ka
    variable_name = 'wse'

variable_to_plot = ds[variable_name]
    fig, ax = plt.subplots(figsize=(15, 10))
contour = ax.contourf(variable_to_plot['x'], variable_to_plot['y'], variable_to_plot, cmap='viridis')
cbar = plt.colorbar(contour)
cbar.set_label('Water Surface Elevation (meters)', fontsize=12)
ax.set_title('SWOT Raster 100m: Lake Mead Reservoir', fontsize=14)

plt.show()

    raster_plot = ds.wse.hvplot.quadmesh(x='x', y='y', rasterize=True, title=f'SWOT Raster 100m: Lake Mead Reservoir')

raster_plot.opts(width=1000, height=600, colorbar=True)

    Creating a Time Series

    SWOT Raster product does not include a time dimension, but it can be inserted by extracting from the file naming convention.

    1. Expand the time range of your earthaccess search to get an adequate range.
    2. Extracting the datetime from the s3 file naming convention then concatenating based on the new time dimension.
    3. Save the output as a new NetCDF file.
    time_results = earthaccess.search_data(
    short_name = 'SWOT_L2_HR_RASTER_2.0',
    bounding_box=(-114.502048,36.060175,-114.390983,36.210182),
    temporal =('2023-12-01 00:00:00', '2024-02-15 23:59:59'),
    granule_name = '*_100m_*',
    count =200
)
    Granules found: 7
    fs_s3 = earthaccess.get_s3fs_session(results=time_results)
    raster = []
for g in time_results:
    for l in earthaccess.results.DataGranule.data_links(g, access='direct'):
        raster.append(l)

print(len(raster))
raster
    7
    ['s3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_007_496_046F_20231211T024534_20231211T024555_PIC0_01.nc',
 's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_008_205_109F_20231221T142037_20231221T142058_PIC0_01.nc',
 's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_008_496_046F_20231231T233040_20231231T233101_PIC0_01.nc',
 's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_009_205_109F_20240111T110542_20240111T110603_PIC0_01.nc',
 's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_009_496_046F_20240121T201546_20240121T201607_PIC0_01.nc',
 's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_010_205_109F_20240201T075048_20240201T075109_PIC0_01.nc',
 's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_010_496_046F_20240211T170050_20240211T170111_PIC0_01.nc']
    def add_time_dimension(ds, file_path):
    # Extract filename from s3 file path
    file_name = file_path.split('/')[-1]
    # Extract date/time string from filename
    date_str = file_name.split('_')[-4][:15]
    # Convert the date string to a datetime object
    time_value = datetime.strptime(date_str, "%Y%m%dT%H%M%S")
    # Assign the time coordinate to the dataset
    ds.coords['time'] = time_value
    return ds
    datasets = []
file_names = []
    for file_path in raster:
    with fs_s3.open(file_path, 'rb') as file:
        file_bytes = file.read()
    file_obj = io.BytesIO(file_bytes)
    dataset = xr.open_dataset(file_obj, engine='h5netcdf')
    dataset_with_time = add_time_dimension(dataset, file_path)
    datasets.append(dataset_with_time)
    file_names.append(file_path.split('/')[-1])
    dataset.close()
    # sorting the time dimension in correct order
datasets.sort(key=lambda ds: ds.time.values)
    result_dataset = xr.concat(datasets, dim='time')
result_dataset
    <xarray.Dataset>
Dimensions:                  (x: 1549, y: 1549, time: 7)
Coordinates:
  * x                        (x) float64 6.788e+05 6.789e+05 ... 8.336e+05
  * y                        (y) float64 3.9e+06 3.901e+06 ... 4.055e+06
  * time                     (time) datetime64[ns] 2023-12-11T02:45:34 ... 20...
Data variables: (12/39)
    crs                      (time) object b'1' b'1' b'1' b'1' b'1' b'1' b'1'
    longitude                (time, y, x) float64 nan nan nan ... nan nan nan
    latitude                 (time, y, x) float64 nan nan nan ... nan nan nan
    wse                      (time, y, x) float32 nan nan nan ... nan nan nan
    wse_qual                 (time, y, x) float32 3.0 3.0 3.0 ... 3.0 3.0 3.0
    wse_qual_bitwise         (time, y, x) float64 8.053e+08 ... 8.053e+08
    ...                       ...
    load_tide_fes            (time, y, x) float32 nan nan nan ... nan nan nan
    load_tide_got            (time, y, x) float32 nan nan nan ... nan nan nan
    pole_tide                (time, y, x) float32 nan nan nan ... nan nan nan
    model_dry_tropo_cor      (time, y, x) float32 nan nan nan ... nan nan nan
    model_wet_tropo_cor      (time, y, x) float32 nan nan nan ... nan nan nan
    iono_cor_gim_ka          (time, y, x) float32 nan nan nan ... nan nan nan
Attributes: (12/49)
    Conventions:                   CF-1.7
    title:                         Level 2 KaRIn High Rate Raster Data Product
    source:                        Ka-band radar interferometer
    history:                       2023-12-15T16:15:45Z : Creation
    platform:                      SWOT
    references:                    V1.1.1
    ...                            ...
    x_min:                         678800.0
    x_max:                         833600.0
    y_min:                         3900500.0
    y_max:                         4055300.0
    institution:                   CNES
    product_version:               01
    ds2 = result_dataset
    variable_name = 'wse'
data = ds2[variable_name]

fig, ax = plt.subplots(figsize=(10, 8))
fig.set_tight_layout({'rect': [0.01, 0.01, 1.0, 1.0]})

contour = ax.contourf(data.isel(time=0), cmap='viridis')
cbar = plt.colorbar(contour)
cbar.set_label('Water Surface Elevation (meters)', fontsize=14) 
times = ds2.time.values

# Function to update the plot for each time step
def update(frame):
    ax.clear()
    contour = ax.contourf(data.isel(time=frame), cmap='viridis',)
    formatted_time = str(times[frame])[:-7]
    ax.set_title(f'Date: {formatted_time}')
    ax.set_xlabel('Longitude', fontsize=14)
    ax.set_ylabel('Latitude', fontsize=14)
    ax.text(0.5, 1.05, 'SWOT Raster 100M Lake Mead Reservoir', transform=ax.transAxes, ha='center', fontsize=14)
    return contour,

# Creating a gif animation
ani = animation.FuncAnimation(fig, update, repeat=True, frames=len(data['time']), blit=True, interval=3000)

output = ('time_series.gif')
ani.save(output, writer='pillow', fps=.5)

with open(output,'rb') as f:
    display(Image(data=f.read(), format='gif'))

plt.close(fig)
ds2.close()
    <IPython.core.display.Image object>

    We can also save the time series to a new netcdf file.

    path = "SWOT_Raster_LakeMead_Time"
os.mkdir(path)
    output_netcdf_path = os.path.join(path, "Output_Time.nc")
result_dataset.to_netcdf(output_netcdf_path)

print(f"Output complete: {output_netcdf_path}")
    Output complete: SWOT_Raster_LakeMead_Time/Output_Time.nc